Our principal question is whether there are fitness differences between HORs, F1s and NORimmigrants. To address this question we will fit a generalized linear model mixed on TLF. In addition to the effect of generation we will also explore covariates that we found were associated with TLF in the McKenzie River reintroduction. These include release length, sex, and julian day of release (fit as a linear continuous variable). We have some priors that suggest we should also explore some interactions: a sex * generation interaction and a sex * length interaction. We also fit year as fixed effect (too few to fit as random like we do in the evaluation: n_year = 4). Finally we will include release group as a random effect. A release group is defined as set of individuals released together at the same location on the same day. Previous analyses have suggested that there is substantial variation in TLF among release groups, and the model fit benefits from the shrinkage associated with including release group.
EDA
Even though we have already done this work for the full dataset, we will conduct a separate exploratory data analysis prior to fitting a model. This will include a careful examination of collinearity, identifying the correct distribution and link function for modeling, and model validation.
Predictors
Let’s explore our predictor variables.
##################################################################
##################################################################
# this is a function to quickly produce biplots between our variables
panel.cor <- function(x, y, digits=1, prefix="", cex.cor = 6)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r1=cor(x,y,use="pairwise.complete.obs")
r <- abs(cor(x, y,use="pairwise.complete.obs"))
txt <- format(c(r1, 0.123456789), digits=digits)[1]
txt <- paste(prefix, txt, sep="")
if(missing(cex.cor)) { cex <- 0.9/strwidth(txt) } else {
cex = cex.cor}
text(0.5, 0.5, txt, cex = cex * r)
}
##################################################################
panel.smooth2=function (x, y, col = par("col"), bg = NA, pch = par("pch"),
cex = 1, col.smooth = "black", span = 2/3, iter = 3, ...)
{
points(x, y, pch = pch, col = col, bg = bg, cex = cex)
ok <- is.finite(x) & is.finite(y)
if (any(ok))
lines(stats::lowess(x[ok], y[ok], f = span, iter = iter),
col = 1, ...)
}
##################################################################
panel.lines2=function (x, y, col = par("col"), bg = NA, pch = par("pch"),
cex = 1, ...)
{
points(x, y, pch = pch, col = col, bg = bg, cex = cex)
ok <- is.finite(x) & is.finite(y)
if (any(ok)){
tmp=lm(y[ok]~x[ok])
abline(tmp)}
}
##################################################################
panel.hist <- function(x, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col="white", ...)
}
F12_mmdata %>%
select(tlf, sex, jday, generation, year, length) %>%
#mutate(tlf = log(tlf)) %>%
pairs(., lower.panel = panel.cor, diag.panel = panel.hist, upper.panel = panel.smooth2)
In the figure above a histogram of each variable is along the diagonal. For the unlabeled factor vairables the levels are (in order) Sex(F, M), Generation (HOR, F1, NORimmigrant).
There are a handful of relationships that jump our even with a simple correlation. The most potentially problematic, as we already noted, is the relationship between generation and year. Generation is also weakly related to release day, sex and length. We will need to explore these in more detail than with this simple figure.
Generation x Release Day
Let’s explore the relationship between release day and generation.
ggplot(data = F12_mmdata)+geom_density(aes(x = jday, fill = generation, color = generation), alpha = 0.4)+theme_bw()+scale_fill_manual(values = c("#228833", "#CCBB44", "#AA3377"))+scale_color_manual(values = c("#228833", "#CCBB44", "#AA3377"))+xlab("Julian Day of Release")
Yes, there is a relationship between generation and release day. Fortunately there is a lot of overlap so our model may be able to parse the effects and collinearity will not prevent us from making clear inferences.
Generation x Sex
ggplot(data = F12_mmdata)+geom_bar(aes(x = generation, fill = sex), position = "fill")+scale_fill_bright()+theme_bw()+ylab("Proportion")
HORs tend to be female biased relative to F1s and NORimmigrants. We will also need to be thoughtful about the relationship between these variables and TLF
Generation x Length
ggplot(data = F12_mmdata)+geom_density(aes(x = length, fill = generation, color = generation), alpha = 0.4)+theme_bw()+scale_fill_manual(values = c("#228833", "#CCBB44", "#AA3377"))+scale_color_manual(values = c("#228833", "#CCBB44", "#AA3377"))+xlab("Length (cm)")
anova(lm(length ~ generation + year, data = F12_mmdata))
plot(emmeans(lm(length ~ generation + year, data = F12_mmdata), "generation", type = "response"))
contrast(emmeans(lm(length ~ generation + year, data = F12_mmdata), "generation", type = "response"), "pairwise")
## contrast estimate SE df t.ratio p.value
## HOR - F1 -2.997 0.367 2721 -8.175 <.0001
## HOR - NORimmigrant -3.715 0.484 2721 -7.679 <.0001
## F1 - NORimmigrant -0.718 0.535 2721 -1.342 0.3722
##
## Results are averaged over the levels of: year
## P value adjustment: tukey method for comparing a family of 3 estimates
HORs are significantly smaller than both F1s and NORimmigrants. F1s and NOR immigrants are the same length.
Distribution
In all recent pedigree work on UWR reintroductions, a negative binomial distribution was a better fit to the data and/or more parsimonious fit than Poisson, QuasiPoisson, and various zero-inflation/hurdle models. However the zero-inflated negative binomial and negative binomial models were pretty close in performance. Let’s re-evaluate given that the zero-inflated model may perform better on this data subset.
For each of the two approaches (negbin and zero-inlfated negbin (zinb)) we’ll conduct model selection then we’ll compare the optimal models. Most of the model selection procedures I already have running won’t work well on the zero-inflated model so we’ll just follow the reccomendations of Zuur et al and rely on backwards stepwise Wald Tests, treating each component (zeros and conditional) separately. We will also ignore the random effects for now.
negbin <- glmmTMB(tlf ~ generation + length + sex+ jday_c+ year, data = F12_mmdata, family = nbinom2)
summary(negbin)
## Family: nbinom2 ( log )
## Formula: tlf ~ generation + length + sex + jday_c + year
## Data: F12_mmdata
##
## AIC BIC logLik deviance df.resid
## 3275.5 3334.6 -1627.8 3255.5 2717
##
##
## Dispersion parameter for nbinom2 family (): 0.733
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.919497 0.591420 -11.700 < 2e-16 ***
## generationF1 0.350368 0.142696 2.455 0.014075 *
## generationNORimmigrant 0.592962 0.152888 3.878 0.000105 ***
## length 0.067292 0.007337 9.171 < 2e-16 ***
## sexM 0.054821 0.093884 0.584 0.559271
## jday_c -0.003360 0.001750 -1.921 0.054776 .
## year2013 0.700092 0.120731 5.799 6.68e-09 ***
## year2014 -0.025644 0.138051 -0.186 0.852637
## year2015 -0.081733 0.154240 -0.530 0.596174
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# drop sex (p = 0.56, wald test)
negbin <- glmmTMB(tlf ~ generation + length + jday_c+ year , data = F12_mmdata, family = nbinom2)
summary(negbin)
## Family: nbinom2 ( log )
## Formula: tlf ~ generation + length + jday_c + year
## Data: F12_mmdata
##
## AIC BIC logLik deviance df.resid
## 3273.8 3327.0 -1627.9 3255.8 2718
##
##
## Dispersion parameter for nbinom2 family (): 0.732
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.874796 0.586633 -11.719 < 2e-16 ***
## generationF1 0.362240 0.141401 2.562 0.0104 *
## generationNORimmigrant 0.607135 0.150995 4.021 5.80e-05 ***
## length 0.067005 0.007324 9.149 < 2e-16 ***
## jday_c -0.003325 0.001750 -1.900 0.0574 .
## year2013 0.697823 0.120719 5.781 7.44e-09 ***
## year2014 -0.027614 0.138022 -0.200 0.8414
## year2015 -0.085282 0.154227 -0.553 0.5803
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# release day is marginal, let's keep it for now
zinb <- glmmTMB(tlf ~ jday_c + sex + generation +length +year , zi = ~ jday_c + sex + generation +length +year , data = F12_mmdata, family = nbinom2)
summary(zinb)
## Family: nbinom2 ( log )
## Formula: tlf ~ jday_c + sex + generation + length + year
## Zero inflation: ~jday_c + sex + generation + length + year
## Data: F12_mmdata
##
## AIC BIC logLik deviance df.resid
## 3260.7 3373.0 -1611.4 3222.7 2708
##
##
## Dispersion parameter for nbinom2 family (): 0.956
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.791119 0.637190 -10.658 < 2e-16 ***
## jday_c 0.001959 0.002184 0.897 0.369537
## sexM 0.052972 0.101544 0.522 0.601902
## generationF1 0.419620 0.149538 2.806 0.005014 **
## generationNORimmigrant 0.522519 0.156951 3.329 0.000871 ***
## length 0.066654 0.007801 8.544 < 2e-16 ***
## year2013 0.969310 0.138809 6.983 2.89e-12 ***
## year2014 -0.058377 0.147297 -0.396 0.691867
## year2015 -0.127272 0.157069 -0.810 0.417771
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.665119 5.431314 -1.411 0.15816
## jday_c 0.139582 0.054149 2.578 0.00994 **
## sexM -0.013771 0.603328 -0.023 0.98179
## generationF1 3.818209 1.431247 2.668 0.00764 **
## generationNORimmigrant 0.234639 1.245773 0.188 0.85060
## length 0.003221 0.044768 0.072 0.94264
## year2013 2.804217 1.530370 1.832 0.06690 .
## year2014 0.891199 1.595979 0.558 0.57657
## year2015 -8.192443 71.075076 -0.115 0.90824
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#sex and length have p-values ~ 1 in the zero part of the model, remove
zinb <- glmmTMB(tlf ~ jday_c + sex + generation +length +year , zi = ~ jday_c + generation +year , data = F12_mmdata, family = nbinom2)
summary(zinb)
## Family: nbinom2 ( log )
## Formula: tlf ~ jday_c + sex + generation + length + year
## Zero inflation: ~jday_c + generation + year
## Data: F12_mmdata
##
## AIC BIC logLik deviance df.resid
## 3256.7 3357.2 -1611.4 3222.7 2710
##
##
## Dispersion parameter for nbinom2 family (): 0.956
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.772952 0.588042 -11.518 < 2e-16 ***
## jday_c 0.001982 0.002099 0.944 0.345063
## sexM 0.054178 0.093352 0.580 0.561666
## generationF1 0.420176 0.147101 2.856 0.004285 **
## generationNORimmigrant 0.522552 0.156919 3.330 0.000868 ***
## length 0.066432 0.007250 9.163 < 2e-16 ***
## year2013 0.968536 0.138544 6.991 2.73e-12 ***
## year2014 -0.059692 0.146460 -0.408 0.683590
## year2015 -0.129179 0.154455 -0.836 0.402958
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.36821 2.40246 -3.067 0.00216 **
## jday_c 0.13934 0.04779 2.916 0.00355 **
## generationF1 3.80170 1.39880 2.718 0.00657 **
## generationNORimmigrant 0.22389 1.22836 0.182 0.85537
## year2013 2.75332 1.29097 2.133 0.03294 *
## year2014 0.84929 1.49261 0.569 0.56936
## year2015 -12.92784 737.07027 -0.018 0.98601
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# next out is sex in the conditional
zinb <- glmmTMB(tlf ~ jday_c + generation +length +year , zi = ~ jday_c + generation +year , data = F12_mmdata, family = nbinom2)
summary(zinb)
## Family: nbinom2 ( log )
## Formula: tlf ~ jday_c + generation + length + year
## Zero inflation: ~jday_c + generation + year
## Data: F12_mmdata
##
## AIC BIC logLik deviance df.resid
## 3255.1 3349.6 -1611.5 3223.1 2711
##
##
## Dispersion parameter for nbinom2 family (): 0.955
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.730949 0.583900 -11.528 < 2e-16 ***
## jday_c 0.002021 0.002104 0.961 0.336654
## generationF1 0.431386 0.145932 2.956 0.003116 **
## generationNORimmigrant 0.536518 0.155205 3.457 0.000547 ***
## length 0.066184 0.007241 9.140 < 2e-16 ***
## year2013 0.966005 0.138511 6.974 3.08e-12 ***
## year2014 -0.063487 0.146531 -0.433 0.664818
## year2015 -0.133281 0.154645 -0.862 0.388769
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.33423 2.40392 -3.051 0.00228 **
## jday_c 0.13963 0.04801 2.908 0.00363 **
## generationF1 3.77084 1.39390 2.705 0.00683 **
## generationNORimmigrant 0.20510 1.22760 0.167 0.86731
## year2013 2.70982 1.28788 2.104 0.03537 *
## year2014 0.78777 1.49189 0.528 0.59747
## year2015 -15.16859 2328.57532 -0.006 0.99480
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# now julian day in the conditional
zinb <- glmmTMB(tlf ~ generation +length +year , zi = ~ jday_c + generation +year , data = F12_mmdata, family = nbinom2)
summary(zinb)
## Family: nbinom2 ( log )
## Formula: tlf ~ generation + length + year
## Zero inflation: ~jday_c + generation + year
## Data: F12_mmdata
##
## AIC BIC logLik deviance df.resid
## 3254.0 3342.6 -1612.0 3224.0 2712
##
##
## Dispersion parameter for nbinom2 family (): 0.922
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.728544 0.583567 -11.530 < 2e-16 ***
## generationF1 0.361337 0.125961 2.869 0.004122 **
## generationNORimmigrant 0.536489 0.153815 3.488 0.000487 ***
## length 0.066156 0.007245 9.131 < 2e-16 ***
## year2013 0.945585 0.136206 6.942 3.86e-12 ***
## year2014 -0.055193 0.144443 -0.382 0.702380
## year2015 -0.143103 0.153281 -0.934 0.350512
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.18147 2.58139 -3.169 0.00153 **
## jday_c 0.14985 0.05019 2.986 0.00283 **
## generationF1 3.98235 1.54102 2.584 0.00976 **
## generationNORimmigrant 0.31799 1.31863 0.241 0.80944
## year2013 3.06659 1.49523 2.051 0.04027 *
## year2014 1.00409 1.73532 0.579 0.56285
## year2015 -15.09148 2346.09962 -0.006 0.99487
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# this looks pretty good
Final ZINB model by backward stepwise selection using Wald Tests is pretty interesting. The conditional part of the model includes an effect of generation, length, and year, while the additional zeros include generation day and year. At face value this suggests that release day affects fitness through the likelihood of reproducing at all, whereas length only matters if you do manage to spawn. This is the basic expectation for how this should work if we assume the zero part of the model largely predicts propensity for pre-spawn mortality whereas the conditional part predicts TLF once you do spawn.
Now let’s compare the two
First we’ll look at model validation
simulateResiduals(negbin, plot = TRUE)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.8833189 0.5585012 0.2131827 0.8928749 0.3074773 0.05568139 0.7841304 0.5851172 0.1772982 0.6566701 0.1105444 0.68838 0.8852255 0.6461571 0.6268861 0.7803337 0.06028185 0.03876886 0.2802067 0.520426 ...
simulateResiduals(zinb, plot = TRUE)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.8860861 0.3415169 0.8554841 0.9076655 0.5378584 0.7534525 0.07024913 0.1368748 0.3072816 0.1536555 0.5930816 0.5762024 0.8107234 0.05952554 0.02173524 0.3002994 0.3933638 0.2059115 0.5031039 0.7490026 ...
Both model fits look good using simulated residuals.
Now let’s compare using AIC, BIC and likelihood ratio tests (the negbin model is nested in teh zinb model)
anova(negbin, zinb)
ZINB best by AIC (delta AIC ~ 20) and likelihood ratio test, but worse by BIC.
So the more complex ZINB model probably provides a better fit to the data. But what do we get for all this added model complexity? Let’s examine the differences in the model fit on the response scale (TLF) using a hanging rootogram.
# refit the final models with a different software to look at rootograms
negbin <- glm.nb(tlf ~ generation +length+ jday_c+year, data = F12_mmdata)
#rescale jday since this software has convergance issues
F12_mmdata %<>% mutate(jday_cs = scale(jday_c))
zinb <- zeroinfl(tlf ~ generation +length +year | jday_cs + generation +year , data = F12_mmdata, dist = "negbin")
rootogram(negbin, main = "Negative Binomial")
rootogram(zinb, main = "Zero-Inflated Negative Binomial")
These are nearly equivalent fits at all TLFs except 1. To me this suggests there isn’t a severe zero-inflation problem and the same processes might drive fitness at both zeros and non-zeros.
Let’s think about it this way, if both models offer similar explanatory power the benefit of the more complex model is additional inference between zero generating and non-zero generation processes (i.e. jday influences PSM whereas length influences TLF once you successfully survive to spawning time). Are we so confident in our model selection procedure that we should stake an entire conclusion on it? How much better is a model with exactly the same variables on both sides of the conditional/zero portions with respect to information/likelihood(AIC, LRT)
zinb <- glmmTMB(tlf ~ generation +length +year , zi = ~ jday_c + generation +year , data = F12_mmdata, family = nbinom2)
# here fit a model where the zero and conditional effects are the same and compare
zinb2 <- glmmTMB(tlf ~ jday_c+ generation +length +year, zi = ~ jday_c + generation +length +year , data = F12_mmdata, family = nbinom2)
# and compare
AIC(zinb, zinb2)
anova(zinb, zinb2)
Delta AIC is ~4, fails to be different in likelihood ratio test. This suggests the more complex model isn’t worth it. Any decent results section would beed to consider both models above (with and without differences between zero and conditional sections).
Instead we should defer to the negbin model. There is nothing wrong with model fit, it is nearly as good of a fit as the zinb, and it is much easier to get our central message across without having to explain the added complexity of model selection in a hurdle/zero-inflation model.
Collinearity
Great, we know we want to use a negative binomial, and we know which variables might be collinear. Let’s see if there is enough information to parse the effects of year, generation and sex on TLF. If we find strong collinearity then our power will be low, our estimates will have large SEs and the inferences we can draw will be limited.
We will do this using the VIF (or more specifically GVIF^(1/2*df)).
# fit beyond optimal model (all possible fixed effect)
beyond_opt <- glm.nb(tlf ~ generation + length + sex+ jday_c+ year , data = F12_mmdata)
vif(beyond_opt)
## GVIF Df GVIF^(1/(2*Df))
## generation 1.988781 2 1.187536
## length 1.215163 1 1.102344
## sex 1.059009 1 1.029082
## jday_c 1.512096 1 1.229673
## year 1.419041 3 1.060065
Very low GVIFs! This is a relief considering the relationship between release day/generation, length/generation and year/generation. Despite these relationships, there is limited correlation between the estimates of their effects and we are likely to have sufficient power/information to identify their independent relationships with TLF.
We are clear to move on to model selection.
Model Selection and Validation
Let’s conduct model selection.
First we will examine the random effect structure, using REML. There is only a single random effect (release group). So, here we compare a mixed model with release group to it’s equivalent glm with no random effects.
beyond_opt_mm <- glmmTMB(tlf ~ generation + length + sex+ jday_c+ year + (1| group ) , data = F12_mmdata, family = nbinom2, REML = TRUE)
beyond_opt <- glmmTMB(tlf ~ generation + length + sex+ jday_c+ year , data = F12_mmdata, family = nbinom2, REML = TRUE)
AIC(beyond_opt_mm, beyond_opt)
BIC(beyond_opt_mm, beyond_opt)
Mixed model is a better performer. We will now switch to model selection for the fixed effect structure and fit using ML. The method is backward stepwise selection using likelihood ratio tests, and a cutoff of p = 0.05.
beyond_opt_mm <- glmmTMB(tlf ~ generation*sex + length+ sex* jday_c+ year + (1| group ) , data = F12_mmdata, family = nbinom2)
drop1(beyond_opt_mm, test = "Chisq")
#drop sex*release day interaction
beyond_opt_mm <- glmmTMB(tlf ~ generation*sex + length + sex+ jday_c+ year + (1| group ) , data = F12_mmdata, family = nbinom2)
drop1(beyond_opt_mm, test = "Chisq")
#drop generation*sex interaction
beyond_opt_mm <- glmmTMB(tlf ~ generation + length + sex+ jday_c+ year + (1| group ) , data = F12_mmdata, family = nbinom2)
drop1(beyond_opt_mm, test = "Chisq")
#drop sex (p value = 0.63)
mm_f12 <- glmmTMB(tlf ~ generation + length + jday_c+ year + (1| group ) , data = F12_mmdata, family = nbinom2)
drop1(mm_f12, test = "Chisq")
# Day of release only marginally improves model fit (delta AIC = 1.8, LRT p value = 0.053). Let's drop it.
mm_f12 <- glmmTMB(tlf ~ generation + length + year + (1| group ) , data = F12_mmdata, family = nbinom2)
drop1(mm_f12, test = "Chisq")
We dropped (in order) sex, and day of release. Now let’s check the fit using simulated residuals.
simulateResiduals(mm_f12, plot = TRUE)
## Object of Class DHARMa with simulated residuals based on 250 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help.
##
## Scaled residual values: 0.03246341 0.4773073 0.1690777 0.7591336 0.1899933 0.3279656 0.7618147 0.4880879 0.1689211 0.5639646 0.07723448 0.6276607 0.9737457 0.679942 0.531472 0.411658 0.8752665 0.815412 0.1754228 0.1984262 ...
Model fit is excellent.
Final Model
summary(mm_f12)
## Family: nbinom2 ( log )
## Formula: tlf ~ generation + length + year + (1 | group)
## Data: F12_mmdata
##
## AIC BIC logLik deviance df.resid
## 3273.1 3326.3 -1627.5 3255.1 2718
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## group (Intercept) 0.03982 0.1995
## Number of obs: 2727, groups: group, 132
##
## Dispersion parameter for nbinom2 family (): 0.763
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.941631 0.590009 -11.765 < 2e-16 ***
## generationF1 0.527212 0.133406 3.952 7.75e-05 ***
## generationNORimmigrant 0.635016 0.163067 3.894 9.85e-05 ***
## length 0.067076 0.007339 9.140 < 2e-16 ***
## year2013 0.707712 0.141223 5.011 5.41e-07 ***
## year2014 -0.027883 0.160516 -0.174 0.862
## year2015 -0.032724 0.183914 -0.178 0.859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Generation (F0, F1, F2 and NORimmigrant), year (fixed effect factor four years) and length significantly improve the fit to the data according to LRT, AIC and Wald tests.
Let’s look at the predicted effects
eff1 <- predictorEffect("generation", mm_f12)
effdf <- as.data.frame(eff1)
effdf$generation <- factor(effdf$generation, levels=c("HOR", "F1", "NORimmigrant")) # relevel the genertions for a nicer plot
#note that this throws an error. w/r/t this error the glmmTMB author (Ben Bolker) states that "the predicted variances are used when computing residuals (which are Pearson residuals by default) for partial residuals plots. I think that if you're not plotting partial residuals, it doesn't matter."
# Since we do not plot partial residuals but instead th95 CI for the predition, we are good here
ggplot(data = effdf, aes(x = (generation), y = fit))+
geom_point(position=position_dodge(width=0.3), size = 3) +
geom_errorbar(aes(ymin = lower, ymax = upper), position=position_dodge(width=0.3), width = 0.1)+ylab("TLF")+xlab("Generation")+theme_bw()
Let’s do some post hoc testing
# we'll use emmeans for this
em <- emmeans(mm_f12, "generation")
contrast(em, "pairwise", adjust = "Tukey", type = "response")
## contrast ratio SE df null t.ratio p.value
## HOR / F1 0.590 0.0787 2718 1 -3.952 0.0002
## HOR / NORimmigrant 0.530 0.0864 2718 1 -3.894 0.0003
## F1 / NORimmigrant 0.898 0.1470 2718 1 -0.659 0.7875
##
## Results are averaged over the levels of: year
## P value adjustment: tukey method for comparing a family of 3 estimates
## Tests are performed on the log scale
The novel finding here is that both NORimmigrants and F1s have greater fitness than HORs, but do not differ from one another.
